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We  liescrlbe  tierr  a  precise  mathematical  fraaevork  for  the  study  of 
discrete-event  systems.  The  idea  is  to  define  a  particular  type  of  sto¬ 
chastic  process,  called  a  generalized  semi-Markov  process,  which  captures 
the  essential  dynamical  structure  of  a  discrete-event  system.  The  paper 
also  attempts  to  give  a  flavor  of  the  qualitative  theory  and  numerical 
algorithms  that  can  be  obtained  as  a  result  of  viewing  discrete-event 
systems  as  GSMP's. 

Keywords:  discrete-event  systems,  simulation,  regenerative  processes,  ' 
Markov  chains,  likelihood  ratios,  variance  reduction,  semi-Markov 
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INTRODDCTIOH 


A  fundamental  obstacle  to  the  study  of  discrete-event  systems  is  the 
lack  of  a  comprehensive  framework  for  the  description  and  analysis  of  such 
systems.  In  this  paper,  we  attempt  to  give  such  a  framework. 

The  idea  that  we  shall  pursue  here  is  to  define  a  particular  type  of 
stochastic  process,  called  a  generalized  semi-Markov  process  (GSMP) ,  which 
captures  the  essential  dynamical  structure  of  a  discrete-event  system.  We 
view  the  GSMP  framework  both  as  a  precise  "language"  for  describing 
discrete-event  systems,  and  as  a  mathematical  setting  within  which  to 
analyze  discrete-event  processes. 

We  start,  in  Section  2,  by  giving  an  abstract  description  of  a 
discrete-event  system.  At  this  level  of  abstraction,  some  of  the  connec¬ 
tions  between  continuous  variable  dynamic  systems  (CVDS's)  and  discrete- 
event  dynamic  systems  (DEBS' s)  become  clear;  this  discussion  is  in  the 
spirit  of  HO  (1987).  Section  3  specializes  the  above  abstract  framework  by 
specifying  a  GSMP  as  a  particular  type  of  event-driven  stochastic  process. 
The  GSMP  structure  is  then  Immediately  applied  to  develop  a  variance  reduc¬ 
tion  technique  that  is  potentially  applicable  to  a  vast  array  of  discrete- 
event  simulations. 

In  Section  4,  the  GSMP  framework  is  specialized  still  further,  thereby 
yielding  the  class  of  time-homogeneous  GSMP's.  These  processes  can  be 
analyzed  via  Markov  chain  techniques.  These  Markov  chain  ideas  are  then 
exploited  in  order  to  obtain  some  qualitative  results  pertaining  to  the 
"long-run"  behavior  of  discrete-event  systems.  Section  5  explores  the 
relationship  between  continuous-time  Markov  chains,  semi-Markov  processes, 
and  GSMP's. 
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Section  6  returns  to  time-homogeneous  GSMP's,  this  time  exploring  the 
qualitative  theory  from  a  regenerative  process  point  of  view.  A  new  type 
of  regenerative  structure  is  described  for  discrete-event  systems  which  are 
scheduled  by  "exponentially  bounded”  distributions.  This  technique  can  be 
applied  to  more  general  GSMP's,  with  some  additional  work,  but  we  present 
here  only  the  current  version.  With  the  aid  of  regenerative  process  ideas  , 
a  strong  law  and  central  limit  theorem  for  discrete-event  systems  are 
established.  In  our  opinion,  these  results  are  typical  of  what  we  can 
expect  to  hold  for  "well-behaved"  discrete-event  processes. 

Finally,  in  Section  7,  we  give  a  flavor  of  the  computational  enhance¬ 
ments  to  discrete-event  simulations  that  are  possible,  by  making  explicit 
use  of  the  GSMP  framework.  Specifically,  likelihood  ratio  Ideas  for 
Importance  sampling  are  briefly  described. 

2.  THE  CVDS/DBDS  ANAL06T 

Suppose  that  we  wish  to  model  the  output  process  (s(t)  :  t  0)  cor¬ 
responding  to  a  (deterministic)  continuous  variable  dynamic  system  (CVDS). 
Frequently,  the  approach  taken  is  to  try  to  represent  s(t)  in  the  form 
s(t)  ®  h(x(t)),  where  x(t)  is  some  suitably  chosen  characterization  of 
the  "internal  state"  of  the  system.  Thus,  given  the  output  process 
(s(t)  :  0  £  8  £  T),  we  can  extend  the  output  process  to  the  Interval 
(T,  T+h)  by  computing  the  internal  state  x  over  the  interval  and  setting 
8(t)  ■  h(x(t)),  T  <  t  £  T  h. 

The  typical  approach  used  to  model  a  discrete-event  dynamic  system 
(DEDS)  is  similar  in  concept.  We  first  recall  that  DEDS  are  frequently 
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used  as  models  of  systems  having  piece-wise  constant  trajectories.  For 
example,  the  trajectory  of  a  queueing  system  Is  constant  between  arrival 
and  departure  epochs  of  customers.  As  a  consequence,  if  (S(t)  :  t  0)  is 
the  output  process  corresponding  to  a  discrete-event  svstem,  it  typically 
takes  the  form 

CO 

(2.1)  S(t)  -Is  I(A(n)  <  t  <  A(n+1))  , 

n-0 

where  we  require  that  0  -  a(0)  <  a(i)  <  ...  (1(A)  represents  an  indica¬ 

tor  function  which  is  one  or  zero  depending  on  whether  or  not  A  occurs.) 

In  the  representation  (2.1),  S  represents  the  output  state  at  the  n'th 

n 

transition  epoch,  and  A(n)  Is  the  instant  at  which  the  n*th  transition 

occurs.  Then,  A  =  A(n+1)  -  A(n)  Is  the  time  between  the  n'th  and 
n+’i 

(n+l)'st  transitions  of  (S(t)  :  t  0). 

To  characterize  the  dynamics  of  the  output  process  (S(t)  ;  t  0) ,  we 

assume  the  existence  of  a  stochastic  sequence  X  -  (X^  :  n  0)  which 

describes  the  time-evolution  of  the  internal  state  of  the  system.  (We 

penult  X  to  be  stochastic  in  order  to  allow  the  discrete-event  system  to 

have  random  behavior.)  We  require  that  the  (S  ,  A  )'8  be  related  to  the 

n  n 

internal  state  sequence  X  via  a  mapping  of  the  form  (S^,  A^)  - 
(h,(X  ),  h«(X  )).  Given  the  output  process  (S(t)  :  0  <  t  <  A(n)),  we  can 
then  extend  the  output  process  to  the  Interval  (A(n),  A(n+1)]  by  first 
computing  We  then  calculate  *  ''2^^n+l^’  A(n+1)  ■  A(n)  + 

and  set  S(t)  -  S(A(n))  for  A(n)  <  t  <  A(n+1).  We  complete  the  extension 
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to  (^(n),  A(rH-l)]  by  calculating  setting  S(A(n+l)) 

=  This  recursive  approach  to  defining  (S(t)  :  t  ^  0)  works, 

provided  that  the  output  process  is  non-explosive  (l.e.,  *(n)  ■*  ®  a.s.  as 
n  ®) . 

The  above  discussion  shows  that  both  the  CVDS  and  DEDS  approaches  to 
modeling  the  output  of  a  system  are,  in  principle,  solved,  once  we 
characterize  the  internal  state  of  the  system.  For  a  CVDS,  perhaps  the 
most  general  characterization  is  to  assume  that  the  Internal  state  process 
X  =  (x(t)  :  t  0)  satisfies  a  relation  of  the  form 

(2.2)  X  =•  f(x) 

for  some  mapping  f.  In  other  words,  for  each  t  0,  this  formulation 
requires  specifying  a  mapping  f(t,.)  for  which 

(2.3)  x(t)  ■  f(t,  x(s)  :  0  £  s  <  ») 

must  hold.  The  analogous  condition  for  a  DEDS  is  to  reaulre  that  there 

exist  a  sequence  of  Independent  r.v.'s  h  “  (h  :  n  >  0)  and  a  map  f 

n  — 

such  that 

(2.4)  X  -  f(X,  h)  . 

This  is  equivalent  to  requiring  the  existence  of  a  family  of  (component) 
mappings  such  that 


(2.5)  .  0  <  k  <  .)  . 

The  resemblance  between  the  equations  (2.2)-(2.3)  and  (2.4)-(2.5)  should  be 
clear. 

Of  course,  the  non-causal  nature  of  (2.3)  creates  difficulties  both 
mathematically  and  computationally.  Furthermore,  formulation  of  a  model 
for  X  directly  in  terms  of  the  mappings  (f(t,*)  :  t  0)  is  often 
unnatural.  As  a  consequence,  it  is  more  typical  to  limit  models  of  the 
internal  state  process  x  *  (x(t)  :  t  0)  to  relations  of  the  form 

x'(t)  “  f(t,  x(s)  :  0  j<  s  <  t) 

(2.6)  s/t 

x(0)  =  Xq  , 

for  some  prescribed  family  of  mappings  (f(t,*)  t  t  >_  0)  and  initial 
condition  Xq.  (Note  that  we  are  now  implicitly  assuming  that  the  Internal 
state  takes  values  in  ]R^.)  The  model  (2.6)  is  causally  defined  in 
terms  of  the  Infinitesimal  characteristics  of  the  system.  The  "local 
specification"  that  is  implicit  in  asserting  that  x  satisfies  a  given 
differential  equation  is  generally  easier  to  formulate  from  a  modeling 
point  of  view  than  the  "global  specification”  implicit  In  (2.3).  Note  that 
the  representation  (2.6)  permits  x  to  be  described  by  differential 
equations  with  delay,  as  well  as  certain  types  of  integro-differential 
equations. 

Of  course,  the  analog  to  the  causal  representation  (2.6)  for  a  DEOS  is 
to  assume  that  the  internal  state  sequence  X  satisfies 

X  *  f  .,(T1^,,  X.  :  0  <  k  <  n) 

(2.7)  ,/t  ^  -  - 
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In  Section  3,  we  describe  a  family  of  discrete-event  systems  In  which  X 
has  the  general  form  given  by  (2.7). 

A  limitation  of  the  causal  models  (2.6)  and  (2.7)  Is  that  the  mathe¬ 
matical  theory  available  to  study  such  unstructured  systems  Is  rather 
poorly  developed.  Foitunately,  In  many  applications  of  CVDS,  It  Is  pos¬ 
sible  to  choose  the  Internal  state  process  x  so  that  the  dynamics  are 
described  by  a  differential  equation  of  the  form 

x'(t)  -  f(t,  x(t)) 

(2.8)  s/t 

x(0)  »  . 


As  we  shall  Indicate  shortly,  the  mathematical  theorv  pertinent  to  (2.8)  is 
quite  extensive. 

The  DEDS  analog  of  the  representation  (2.8)  Is  to  require  that  the 
Internal  state  sequence  X  satisfy  a  recursion  of  the  form 


(2.9) 


s/t 


^iH-1  *  ’^n+l^ 


Xp  . 


Such  representations  can  often  be  obtained  for  systems  satisfying  (2.7), 
provided  that  a  judicious  choice  of  state  space  Is  made.  Specifically,  it 
is  often  possible  to  adjoin  "supplementary  variables"  to  a  state 

descriptor  X^  satisfying  (2.7),  to  obtain  a  new  state  sequence 
-  (X^,  X^^)  satisfying  (2.9). 

The  mathematical  power  of  the  representation  (2.9)  is  a  consequence  of 
the  following  easily  proved  result. 
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PROPOSITION  1.  Suppose  that  X  satisfies  (2.9),  and  that  the  r.v.'s 


{xq,  :  n  1}  are  independent.  Then,  X  is  a  Markov  chain  (i.e., 

P{X  E  •  Ix^,  ...,  X  }  =  P{X  . ,  e  •  |X  }). 
n+1  '0  n  n+1  '  n 


A  substantial  literature  on  the  theory  of  Markov  chains  can  be  applied 
to  the  analysis  of  DEDS  for  which  the  internal  state  sequence  satisfies  the 
conditions  of  Proposition  1.  Similarly,  the  vast  mathematical  theory  on 
differential  equations  is  directly  relevant  to  CVDS  satisfying  the  repre¬ 
sentation  (2.8).  For  example,  existence  theory  for  the  differential 
equation  (2.8)  basically  yields  conditions  under  which  there  exists  an 
output  process  (compatible  with  (2.8))  which  can  be  defined  over  the  entire 
seml-lnf Inlte  Interval  [0,*).  The  DEDS  counterpart  Involves  deriving 
conditions  under  which  (S(t)  ;  t  0)  Is  non-explosive. 

Much  of  the  differential  equations  literature  on  systems  satisfying 
(2.8)  pertains  to  systems  obeying  the  stronger  condition 


(2.10) 


s/t 


x’(t) 

x(0) 


f(x(t)) 


This  literat'j'^?  typically  focuses  on  the  large-time  behavior  of  the  inter¬ 
nal  state  process  (x(t)  ;  t  0).  This,  in  turn,  is  strongly  related  to 

the  study  of  the  set  (x  :  f(x)  •  0}  of  equilibrium  points  for  (2.10). 

The  DEDS  counterpart  to  (2.10)  requires  a  model  formulation  in  which 
the  Internal  state  sequence  X  takes  the  form 


(2.11) 


s/t 


^n+1  *  ^n+1^ 


‘o  • 


The  following  result  is  easily  demonstrated,  and  so  the  proof  is  omitted. 
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PIDPOSITION  2.  Suppose  that  X  satisfies  (2.11).  In  addition,  assume 
that  ;  n  l)  Is  a  collection  of  independent  identically  distributed 

(i.i.d.)  r.v.'s,  independent  of  Xq.  Then,  X  is  a  time-homogeneous  Markov 
chain  (l.e.,  there  exists  a  transition  function  P(x,  A)  such  that 

"■'Vl  '  •  1*0 . *n>  ■ 

As  in  the  CVDS  setting,  much  of  the  mathematical  literature  on  Markov 
chains  of  the  form  (2.11)  concentrates  on  study  of  the  long-run  behavior  of 
the  system.  The  concept  of  equilibrium  point  is  now  replaced  by  that  of  an 
invariant  probability  distribution.  A  probability  distribution  x  is  said 
to  be  invariant  for  the  (time-homogeneous)  Markov  chain  X  if 

(2.12)  x(dy)  -  J  x(dx)  P(x,  dy) 

Z 

(Z  is  the  state  space  of  X).  In  the  presence  of  irreduclbllity  hypothe¬ 
ses  on  X,  the  existence  of  an  invariant  probability  distribution  x 
typically  Implies  that  for  each  (measurable)  subset  A  of  Z, 

^  I  KX.  e  A)  -*■  x(A)  a. 8. 

”  k-0 

as  n  ®,  for  every  possible  initial  condition  Xp.  The  analogous  CVDS 
concept  is  that  of  (global)  stability:  there  exists  x  such  that 
x(t)  X  as  t  •♦■  *,  for  every  possible  initial  condition  Xp.  See 
CODDINGTON  and  LEVINSON  (1955)  for  further  details. 

As  the  above  discussion  suggests,  most  of  the  mathematical  theory 
pertaining  to  systems  of  the  form  (2.10)  and  (2.11)  is  qualitative  in 
nature.  A  major  computational  difference  between  CVDS  and  DEDS  is  that  the 
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numerical  determination  of  equilibrium  points  Is  significantly  simpler  than 
that  of  calculating  Invariant  probability  distributions.  Nevertheless,  It 
Is  our  view  that  the  above  analogy  between  CVDS  and  DEDS  Is  useful  In 
developing  -n  understanding  of  the  major  theoretical  Issues  arising  from 
discrete-event  systems. 


3.  CKMERALIZED  SEKl-MARKOV  PROCESSES 

Consider  a  discrete-event  system  In  which  the  Internal  state  sequence 
X  has  the  causal  representation  (2.7).  Then,  If  ***•  ^n^ 

the  o-fleld  corresponding  to  the  history  of  X  up  to  time  n,  we  find 
that 


(3.1) 


PiX  . ,  e 
n+1 


!  -  P( 


^n> 


where  the  conditional  probability  appearing  on  the  right-hand  side  of  (3.1) 
Is  defined  by 


(3.2)  P(*  ;  Xp,  ...,  x^)  -  ’^0’  V  * 

The  discrete-event  system  evolves  In  time  by  recursively  generating  X  . , 

n+ 1 

from  the  conditional  distribution  (3.2).  Once  is  obtained, 

and  can  be  calculated  by  using  the  transformations  hj^  and  h2. 

Of  course,  most  discrete-event  systems  of  Interest  have  more  structure 
than  that  which  we've  described  above.  Specifically,  discrete-event 
systems  are  typically  characterized  by  two  different  types  of  entitles, 
namely  states  and  events.  In  a  queueing  network.,  the  states  generally 
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correspond  to  queue-length  vectors  describing  the  number  of  customers  at 
each  station  of  the  network.  The  set  of  events  lists  the  different  ways  In 
which  the  queue-length  vector  can  change  as  a  customer  completes  service  at 
one  station  and  moves  to  the  next. 

To  mathematically  describe  the  dynamics  of  this  type  of  system,  we  let 
S  denote  the  (finite  or  countable)  set  of  states  and  E  denote  the 
(finite  or  countable)  set  of  events.  A  state  s  e  S  is  termed  a  physical 
state ,  In  order  to  distinguish  such  states  from  the  state  space  correspond¬ 
ing  to  the  internal  state  sequence  of  the  discrete-event  system.  For  each 
8  e  S,  let  E(s)  be  a  non-empty  finite  subset  of  E  denoting  the  set  of 
events  that  can  trigger  transitions  out  of  state  s. 

EXAMPLE  1:  To  model  an  open  queueing  network  with  d  stations  and  one 
class  of  customers,  we  let  S  •  •••  ^  (d  times).  The  vector 

s  •  (s(l),  ...,  s(d))  e  S  will  then  represent  the  queue-lengths  (Including 
the  customer  at  the  server)  at  each  of  the  d  stations.  A  state  transi¬ 
tion  occurs  via  either  of  the  following  possibilities:  an  external  arrival 
event  or  a  departure  event.  Thus,  E(s)  *  {(1,  1)  :  1  £  i 

J  {(1,  2)  :  1  £  1  £  d,  s(l)  _>  1} ,  where  (1,  1)  corresponds  to  an  external 
arrival  to  station  1  and  ( j ,  2)  denotes  a  departure  from  station  j. 

For  each  event  e  t  E(8),  we  can  associate  a  clock.  The  reading  c^ 
on  clock  e  can  be  viewed  as  representing  (In  some  rough  sense)  the  amount 
of  time  that  has  passed  since  clock  e  was  last  activated. 
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EXAMPLE  1  (continued):  Suppose  that  there  exists  a  single  server  at  each 
of  the  d  stations;  each  server  works  at  unit  rate.  For  e  «  (1,  1),  c^ 
corresponds  to  the  amount  of  time  that  has  elapsed  since  the  last  external 
arrival  to  station  1.  For  e  ( j ,  2),  c^  denotes  the  amount  of  time 

that  has  passed  since  service  was  Initiated  on  the  customer  at  the  server 
at  station  J. 

The  clock  readings  for  event  e  Increase  at  a  speed  Thus,  the 

rate  at  which  c^  Increases  may  depend  on  the  physical  state  occupied  by 
the  discrete-event  system. 


EXAMPLE  1  (continued):  If  each  of  the  single  servers  serves  at  unit 
rate,  then  r  »  1  for  all  s  e  S,  e  e  E(s),  On  the  other  hand,  to  model 
a  server  for  which  the  service  rate  Is  proportional  to  the  number  of  cus¬ 


tomers  in  queue  at  the  station,  we  would  set  r  ■  r,  •  s(j)  for 

SG  j 

e  =  (j,  2). 


We  let  =  {-1}  ‘J  (0,  “)  and  assume  that  c^  t  We 

adopt  the  convention  that  If  c^  «  -1,  then  e  Is  currently  not  active. 

(For  example,  event  ( j ,  2)  Is  not  active  when  s(j)  ■  0.)  Thus, 

+  E 

C(s)  =  {c  e  (R  ,)  :  c  ■  -1  Iff  e  ^  E(s)}  Is  the  set  of  clock 

-1  e 

readings  possible  In  s  e  S. 

We  now  give  a  rough  description  of  the  dynamics  of  the  discrete-event 
system.  Suppose  that  at  time  ^(n),  the  system  has  just  entered  state 
8  e  S  and  that  the  clock  reading  at  that  Instant  is  represented  by 
c  e  CCs).  Each  clock  e  e  Eis)  will  now  compete  to  trigger  a  transition 
out  state  s.  The  system  evolves  by  first  probabilistically  generating. 
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for  each  event  e  e  E(s),  a  “residual  lifetime"  r.v.  which  represents  the 

amount  of  time  remaining  until  event  e  triggers  a  transition  out  of  state 

8.  This  residual  lifetime  has  a  distribution  which  depends,  of  course,  on 

the  values  of  c^  and  (It  may  also  depend  on  the  history  ^  in  a 

more  complicated  way.)  Having  generated  residual  lifetimes  for  all  the 

events  e  e  E(s),  the  trigger  event  e*  is  simply  the  next  event  to  be 

scheduled  by  the  DEDS.  In  other  words,  the  trigger  event  is  the  event  e 

having  the  minimal  residual  lifetime.  This  minimal  residual  lifetime 

therefore  yields  A  (l.e.,  the  time  between  the  n'th  and  (n+D'st 

n+1 

transitions) . 


A  new  physical  state  s*  is  now  chosen  stochastically;  its  distribu¬ 
tion  typically  depends  on  both  the  previous  state  s  and  the  trigger  event 
e*.  The  events  e  e  0(s';  s;  e*)  =  E(s')  -  (E(s)  -  {e*})  have  their  clock 
readings  incremented  appropriately  to  refect  the  speed  r^^  and  the 
passage  of  time  (Clocks  e  e  0(s';  s,  e*)  are  "old”  clocks  that 

continue  to  run  in  state  s'.)  Clearly,  the  events  e  e  N(s';  s,  e*) 

=  E(s')  -  0(s';  s,  e*  are  the  "new”  clocks,  which  will  satisfy  c  ,  “0 

8  6 

at  time  A(n+1).  Thus,  we  have  calculated  the  physical  state  and  clock 
readings  at  time  ^(n+l).  The  process  can  now  be  repeated  recursively  to 
obtain  the  physical  state/clock  readings  at  A(n+2),  A(n+3),  ...  . 

We  shall  now  describe  the  time  evolution  of  the  system  more  pre¬ 
cisely.  Let  Z  -  'J  {s}  [0,  ®)  C(s).  The  Internal  state  sequence  X 

stS 

will  be  assumed  to  take  values  in  Z.  Specifically,  the  internal  state  at 


transition  n  is  given  by  X  »  (S  ,  A  ,  C  )  e  Z.  Thus,  the  first  two 

n  n  n  n 
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components  of  correspond  to  the  physical  state  and  "holding  time"  of 

the  system  at  the  n'th  transition.  The  third  component  Is  a  vector 

denoting  the  state  of  the  clocks  at  time  ^^Cn).  Note  that  the  mappings 

hj,  h^  of  Section  2  are  given  by  hj(8,  t,  c)  ■  s,  h2(8,  t,  c)  ■  t. 

If  e  E,  1  0,  the  vector  ■  (xq . x^)  Is  a  possible 

realization  of  X  ■  (X„,  X  ).  We  now  need  to  describe  the  condl- 

n  u  n 

tional  probability  (3.2)  (i.e.,  p{X  ,,e  •  |X  ■  x  })  In  more  detail.  To 

n'+’i  n  n 

set  the  stage,  we  assume  that  for  all  s  c  S,  there  exists  an  event 

e  e  ECs)  such  that  r  >0.  Thus,  In  every  state,  there  exists  at  least 

one  clock  with  positive  speed.  We  further  let  p(s';  x  ,  e)  be  the 

n 

(conditional)  probability  that  S  . ,  equals  s' ,  given  that  X  equals 

n+t  n 

x^  and  the  trigger  event  transition  n  +  1  is  e.  Also,  we 

assume  the  existence  of  a  family  F(*;  x^,  e)  of  probability  distribu¬ 
tion  functions  such  that  F(0;  x^,  e)  •  0.  (i.e.,  F(*;  x^,  e)  corre¬ 

sponds  to  a  positive  r.v.)  The  distribution  ?(•;  x^,  e)  helps  govern 
the  "residual  life"  of  clock  e,  given  that  X^  ■  x^.  We  require  that 

for  each  x  ,  there  exists  at  most  one  clock  e  such  that  F(t;  x  ,  e) 
n  n 

is  not  continuous  as  a  function  of  t.  This  guarantees  that  the  trigger 
event  e*  ,  will  be  uniquely  defined  for  each  x  .  Set 


F^(x;  x^,  e) 


F(ax;  x^,  e) ,  a  ^  0  , 


F,(*;  e) 


1  -  F  (x;  X  ,  e) 
a  n 


G(dt;  x^,  e) 


-  F  (dt;  X  ,  e)  IT  F  (t;  x  ,  e') 

8  ,e  e'eE(s  )  s  ,e* 

n’  n  n’ 

e'  ^  e 
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(s  is  the  first  component  of  x  ).  Note  that  G(dt;  x  ,  e)  represents 
n  n  n 

the  probability,  conditional  on  x  ,  that  A  e  dt  and  e*  ,  *  e.  We 

n  n+i  n+i 

can  now  rigorously  define  the  conditional  probability  structure  of  the 

Internal  state  sequence  X.  For  A  -  {s' }  [0,  T]  X  (-®,  a  ],  let 

ecE  ® 


(3.3) 


■* 


P{X_^,  e  A  I  X  -  X  } 
n+1  '  n  n 


I  P(8’;  e*) 


e*eE(s  ) 
n 


n  I(a  ,  >  -1) 
e’^ECs')  ® 


n  I(a  ,  0)  •  /  G(dt;  X  ,  e*) 

e'eN(s' ;8^,e*)  ®  [0,T)  ” 


n 

e'e0(s' ;8^,e*) 


+  t  •  r^  ) 


e  —  e'  ,n 


.e- 


(c^  is  the  clock  reading  vector  of  x^).  The  product  of  indicators  over 
e'  i  ECs'Xe'  e  N(s';  s^,  e*))  represents  the  fact  that  clocks 
e  ^  E(s')(e'  e  NCs';  s^,  e*))  have  clock  readings  of  -1(0).  The 
product  of  indicator  functions  over  0(s';  s^,  e*)  corresponds  to  the  fact 
that  the  "old"  clocks  need  to  be  properly  Incremented  to  their  new  values 
at  A(n+1). 

The  conditional  probability  distribution  (3.3)  asserts  that  we  may 

generate  X  .  ,  from  X  in  the  following  way.  We  first  generate  Inde- 
n+l  n 

pendent  r.v.'s  Y  (eeE(S))  from  the  conditional  distributions 
e  ,n  n 

F(*;  X  ,  e).  Then,  A  ■  mln{Y  /r„  ;  e  t  E(S  )}  and  the  trigger 

n  nri  e  ,n  a  ,e  n 

’  n 

event  e*  ,  is  the  (unique)  event  e  £  E(S  )  which  achieves  the  minimum 

ttr  1  n 

for  We  then  generate  from  p(»;  X^,  Finally,  we  set 
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-1 


^n+l,e 

^n+l,e 


C 

n,e 


+  A 


n+1 


’  *  '  ®n’  ®otl'  ■ 


'S  "  “<Vl'  *n-  'ItP  • 

•%  » 


We  call  a  discrete-event  system  with  an  Internal  state  sequence  X 

satisfying  (3.3)  a  ( t Ime-lnhomogeneous )  generalized  semi-Harkov  process. 

The  term  "time-lnhomogenelty"  reflects  the  fact  that  the  "residual  life" 

distributions  and  state  transition  probabilities  p(s';  x  ,  e)  can 

n 

depend  explicitly  on  the  entire  history  of  X.  For  example,  these  proba¬ 
bility  distributions  may  depend  explicitly  on  A(n)  (i.e.,  the  time  of  the 
n'th  transition);  see  Section  7  for  further  details.  Furthermore,  as  we 
shall  show  in  Section  5,  these  processes  do  indeed  extend  the  notion  of 
semi-Markov  process,  thereby  justifying  use  of  the  term  generalized 
semi-Markov  process.  For  the  remainder  of  this  paper,  we  will  refer  to 
discrete-event  systems  satisfying  (3.3)  simply  as  generalized  semi-Markov 
processes. 

We  conclude  this  section  with  an  illustration  of  how  we  can  exploit 
the  causal  structure  of  DEDS  satisfying  (2.7)  to  obtain  improved  statisti¬ 
cal  efficiency  for  associated  simulations.  Specifically,  suppose  that  we 
wish  to  calculate,  via  simulation,  an  expectation  of  the  form 


a  -  Ef(XQ, 
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(f  :  Z  •••  ^  Z  (n+1  times)  TR).  The  standard  approach  would  first 

replicate  the  r.v.  f(X„,  X  )  m  independent  times,  and  then  form  the 

u  n 

sample  mean  of  the  m  observations. 

However,  an  alternative  estimator,  based  on  control  variates,  is  often 
available.  Suppose  there  exists  a  random  d-vector  X  such  that  Ex  ■  0. 

Such  a  vector  X  is  known,  in  the  simulation  literature,  as  a  control.  Then, 

C(X)  -  f(XQ,  ...»  X^)  -  X*^x 

is  an  unbiased  estimator  of  a  for  all  X  e  (We  adopt  the  convention 

that  all  elements  of  are  written  as  column  vectors) .  Since  X  is 

at  our  disposal,  we  may  choose  X  to  minimize  var  C(X).  The  optimal 
value  of  X  is  given  by 

(3.4)  X*  -  (EXX^)“^  •  EXfCXp . X^)  . 

To  implement  the  method  of  control  variates,  we  generate  m  independent 

copies  of  the  pair  (f(X,,,  ...,  X  ),  x)»  If  X  is  a  sample-based  esti- 

u  n  m 

mate  for  X*,  we  obtain  an  (asymptotic)  improvement  over  the  original  esti¬ 
mator  by  using  a  sample  mean  of  the  C(X  )'8  rather  than  f(X-,  ...,  X  )’s. 

m  U  n 

The  basic  idea  underlying  the  use  of  control  variates  is  that  we  are 
"filtering  out"  the  noise  in  ^(Eq,  •••,  X^)  due  to  X;  this  then  reduces 
the  variance  of  the  resulting  estimator. 

The  key  to  the  method  of  control  variates  is  to  obtain  an  easily 
calculated  control  X  which  is  highly  correlated  with  ^(Xq,  ...,  X^).  It 
turns  out  that  the  causal  structure  represented  by  (2.7)  can  be  used  to 
easily  obtain  control  variates.  Suppose  that  the  conditional  distribution 
(2.7)  has  the  property  that  for  some  real-valued  function  g,  the  condi¬ 
tional  mean 


16 


may  be  easily  calculated.  For  example.  If  the  DEDS  Is  a  GSMP,  the  condi¬ 
tional  means  e  ^^n  *  ^n^  often  have  simple  analytical  closed- 

forms. 

Let  -  g(X^)  -  g(X  ,)  for  ra  >  1.  If  Eg^(X  )  <  ®  for  n  >  0, 
n  n  n—  l  —  n  — 

it  may  be  easily  verified  that  the  r.v.'s  D, ,  . . D  are  martingale 

1  z  n 

differences  with  respect  to  the  sequence  of  o-fields  (  .^  :  m  >  0).  Con- 

IB  ■“ 

sequently,  ...,  are  orthogonal  mean-zero  r.v.'s.  Since  the 
have  mean  zero,  it  follows  that  X  “  (D^,  ...,  Is  a  control.  The 

orthogonality  of  the  D^'s  Implies  that  EXX^  is  a  diagonal  matrix.  Thus, 

(3.4)  simplifies  to 

(3.5)  X*  -  ED^f(XQ,  ...,  X^)/ED^  , 

2  2 

provided  that  Ef  (X^,  ...,  X^)  <  ®  and  ED^^  >  0.  Hence,  an  advantage 

of  the  martingale  controls  described  here  is  that  X  need  estimate  only 

2 

2n  parameters  in  (3.5),  as  opposed  to  (n  +  3n)/2  in  (3.4), 

The  above  discussion  shows  that  the  method  of  control  variates  is 
generally  applicable  to  DEDS  in  which  the  internal  state  sequence  is 
causally  generated.  In  particular,  martingale  control  variate  schemes  can 
be  applied  to  GSMF's. 
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4.  TIME-HOMOGENBOUS  GSMP'S 

In  this  section,  we  examine  a  class  of  GSMP’s,  which  also  satisfy 
(2.11).  As  shown  in  Proposition  2,  this  will  guarantee  that  X  is  a  time- 
homogeneous  Markov  chain. 

In  many  discrete-event  systems,  the  constituent  conditional  probabil¬ 
ity  distributions  F(*;  x^,  e)  and  p(»;  x^,  e)  defining  a  GSMP 
simplify  considerably. 

EXAMPLE  1  (continued).  Suppose  that  customers  are  routed  through  the 
queueing  network  via  a  substochastic  Markovian  switching  matrix  P.  Then, 
the  state  transition  probabilities  p(*;  x^,  e)  take  the  form 
p( • ;  s^,  e) .  This  means  that  the  probability  distribution  of  depends 

only  on  the  current  physical  state  s^  and  the  trigger  event  e.  Let  e^ 
denote  the  i'th  unit  vector.  The  specific  form  of  the  state  transition 
probabilities  is  given  by 

1  if  s'  -  8  +  e. 

p(s';  8,  (i,l))  - 

0  if  s'  ^  8  +  e^ 

P. .  if  s'  -  s  +  e.  -  e.  ,  8(1)  1 

p(8';  s,  (1,2))  -  ^  ^ 

d 

M  -  I  p..  if  s'  -  8  -  e^,  8(1)  2  1  . 

J»1  ^ 

Suppose  that  we  further  assume  that  the  external  arrival  stream  to  the 
i'th  station  is  a  renewal  process  with  continuous  inter-arrival  distribu¬ 
tion  F^.  Also,  suppose  that  each  server  employs  a  first-come/ first-serve 
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queueing  discipline  in  which  the  service  requirements  for  the  consecutive 
customers  served  at  the  i'th  station  are  1.1. d.  with  common  continuous 
d is  t  r ibut ion  . 

If  e  ■  (i,l),  the  clock  reading  c  is  the  amount  of  time  (at  the 

instant  ^(n))  that  has  passed  since  the  last  customer  arrived  externally 

to  station  i.  For  e  ■  (j,2),  the  clock  reading  c  is  interpreted  as 

n,e 

the  amount  of  service  requirement  that  has  been  processed  on  the  customer 
that  is  in  service  at  station  j  at  time  A(n). 

Assume  that  the  service  rate  of  the  server  at  station  i  is 
r^^  •  s(i),  so  that  the  rate  is  proportional  to  the  number  of  customers  in 
queue  at  station  i.  Given  the  above  assumptions,  the  conditional  distri¬ 
bution  P(dt;  X  ,  e)  “  F(dt;  c  ,  e)  so  that  the  conditional  proba- 
n  n,e 

blllty  distribution  function  for  clock  e  depends  on  the  history  of  X 
only  through  c^  For  e  -  (1,1),  the  exact  form  of  the  conditional 
distribution  is  given  by 


F(t;  c 


n.e 


e) 


F,(t 


c  -)/F.(c 
n  ^6  X  R  j 


whereas  for  e  ■  (j,2). 


Building  on  the  above  example,  suppose  that  we  have  a  GSMP  for  which 
there  exists  a  family  of  distributions  (F^  :  e  c  E)  and  a  family  of 
state  transition  probabilities  (p(*;  s,  e)  ;  s  e  S,  e  e  E(s))  such  that 
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(4.1) 
and 

(4.2) 


p(*;  8^,  e) 


If  a  GSMP  satisfies  the  additional  conditions  (4.1)-(4.2),  ue  refer  to 

the  discrete-event  system  as  a  time-honogeneous  GSMP.  Noting  that 

p{X^,  e*|X  -x}  can  be  represented  in  the  form  P(x  ,  •),  we  see 
n+1  '  n  n  n 

that  the  internal  state  sequence  of  a  time-homogeneous  GSMP  is  a  time- 
homogeneous  Harkov  chain. 

The  Markov  structure  of  a  time-homogeneous  GSMP  can  be  fruitfully 
exploited  to  study  the  long-run  behavior  of  the  corresponding  discrete- 
event  systeci.  The  following  result  shows  that  the  long-run  behavior  of  the 
output  process  (S(t)  :  t  0)  can  be  calculated  from  that  of  the  internal 
state  sequence  X.  (In  order  that  Theorem  I  hold,  the  system  need  not  even 
be  a  GSMP.) 


THEOREM  1.  Let  f  be  a  real-valued  function.  Suppose  that 


n-1 


1  *  a.s. 


j-0 


y  j+1 


1  a.s. 


where 


n  -> 


^1*  ^2’  finite  r.v.'s  with  ^2  ^  0  a.s.  Then, 
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t 


I  f(S(u))du 
0 


2 


as  t  ->•  ®. 


The  proof  of  this  result  Is  similar  to  that  of  Proposition  2  of  GLTNN 
and  IGLEHART  (1988a)  and  so  is  omitted. 

To  examine  the  long-run  behavior  of  the  Harkov  chain  X,  we  study  the 
question  of  existence  of  Invariant  probability  distributions  for  X. 

TBBOSEM  2*  Let  X  be  the  Markov  chain  corresponding  to  a  time- 
homogeneous  GSMP  with  the  following  properties: 


1) 

|S|  <  » 

11) 

F  Is  continuous 
e 

with  F^(0)  ■  0,  for 
e 

all  e  e  E 

ill) 

r  >  0 
se 

for  all 

8  e  S,  e  e  E(8) 

Iv) 

F^(c)  < 
e 

1  for  all 

e  e  E,  c  <  ® 

v) 

For  all 

e  >  0,  e 

e  E,  there  exists  R 

such  that 

F  (K+c)/F  (c)  <  e 

9  0 

uniformly  In  c. 

Then,  X  has  an  Invariant  probability  distribution  x. 

For  the  proof,  see  the  appendix.  Suppose  that  F^(*)  has  a  density 

f  (*).  Then  Che  hazard  rate  function  h.(*)  Is  defined  via  h  (t) 
e  .  ' ■  ■ '  -  '  e  e 

■  f  (t)/F  (t)  (F  (t)  •  1  -  F  (t)).  Condition  v)  Is  satisfied  if 
e  e  e  e 

the  hazard  rate  function  Is  bounded  below  by  a  positive  constant  (l.e., 
inf{h^(t)  :  t  0}  >  0).  It  Is  also  satisfied  by  any  finite-mean  dis¬ 
tribution  which  Is  new  better  than  used  in  expectation  (see  BARLOW  and 
PROSCHAN  (1975)  for  a  discussion  of  such  distributions). 
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Hypotheses  i)  and  v)  of  Theorem  2  guarantee  that  the  Markov  chain  X 
spends  a  large  fraction  of  time  in  compact  subsets  of  I;  this,  in  turn, 
guarantees  the  type  of  "positive  recurrence”  needed  to  obtain  the  existence 
of  invariant  probability  measures.  The  proof  also  demands  that  the  transi¬ 
tion  function  of  X  be  continuous  in  a  certain  sense;  conditions  ii)-iv) 
yield  the  required  continuity. 

Let  P  (•)(£(•))  denote  the  probability  (expectation)  on  the  path 

I*  r* 

space  of  X  under  which  X^  has  distribution  p. 

THEOREM  3.  Let  X  be  the  Markov  chain  corresponding  to  a  t ime- homogene¬ 
ous  GSMP  with  the  following  properties: 

i)  |S(  <  «, 

ii)  /  t  F  (dt)  <  «,  e  e  E, 

(0,«)  ® 

iii)  For  all  s,  s'  e  S,  E  e  ECs)  with  pCs’;  s,  e)  >  0,  there 

exists  e'  e  N(s';  s,  e)  such  that  r  ,  ,  >  0. 

s  ,e 

Then,  if  X  has  an  invariant  probability  distribution  n,  E^(f(sp)iAj  <  ®. 

For  the  proof  of  Theorem  3,  see  the  appendix.  The  point  of  Theorem  3 
is  that  it  gives  sufficient  conditions  for  the  finiteness  of  E„A,  and 
E^|f(sQ)(A^.  Such  moment  conditions  are  necessary  in  order  to  apply  the 
ergodic  theorem.  The  following  result  is  an  immediate  consequence  of 
Birkhoff's  ergodic  theorem  and  Theorem  1.  (We  need  the  continuity  of  the 
F  (•)'s  in  order  to  guarantee  that  P{A  >  0(x^  ■  x)  ■  1  for  all  x  €  E. 
This  ensures  that  >  0  a.s.). 
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PROPOSITION  3.  Let  X  be  the  Markov  chain  corresponding  to  a  time- 

homogeneous  GSMP.  Suppose  that  there  exists  an  invariant  probability 

distribution  x  for  X  such  that  A,  <  ®  and  E  IfCS^)  lA,  <  We 

X  1  X '  0  1 

further  assume  that  F  (•)  is  continuous  for  all  e  e  E.  Then, 

e  ’ 


■p  /  f(S(u))  du 

0 


as  t  ->■  ®,  where  ^  is  the  invariant  ff-field  corresponding  to  X. 

Theorems  1  to  3,  together  with  Proposition  3,  give  conditions  under 
which  a  discrete-event  system  "settles  down"  to  a  steady-state.  It  should 
be  emphasized  that  the  results  merely  assert  existence  of  a  steady-state 
and  say  nothing  about  uniqueness.  In  particular,  under  the  conditions 
given  above,  it  is  quite  possible  for  the  system  to  have  multiple 
steady-state  distributions.  The  particular  steady-state  distribution 
governing  the  discrete-event  system  then  depends  on  the  initial  state  X^. 

Related  results  on  existence  of  invariant  probability  measures  for 

time-homogeneous  GSMP's  appear  In  KONIG,  MATTHES,  and  NAWROTZKI  (1967), 

(1974)  and  WHITT  (1980).  The  latter  paper  also  gives  conditions  under 

which  the  invariant  probability  distribution  x  is  continuous  in  the 

state-transition  probabilities  p(s';  s,  e)  and  distributions 

(F^  :  e  e  E). 
e 

In  Section  6,  we  return  to  this  steady-state  theme.  The  regenerative 
machinery  used  there  establishes  both  existence  and  uniqueness  results  (but 
under  different  conditions  than  those  discussed  here). 
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5.  CCWTIIIDOUS-TIME  MASKOV  CHAINS  AND  SEKE-NARKOV  PROCESSES 

Our  objective  here  Is  to  briefly  Indicate  some  of  the  connections 
between  continuous-time  Markov  chains,  semi-Markov  processes,  and  GSMP's. 

Basically,  any  time- homogeneous  GSMP  in  which  all  the  event  distribu¬ 
tions  Fg(*)  are  exponential  is  a  continuous  time  Markov  chain.  More 
precisely,  consider  a  time-homogeneous  GSMP  for  which 

F^(dt)  -  X(e)  exp(-\(e)t)dt 

(^.(e)  >  0)  for  all  e  e  E.  Then,  the  conditional  probability  distribu¬ 
tions  defined  by  (4.2)  take  the  form 

(5.1)  F(t;  X  ,  e)  =  exp(-X.(e)t) 

n 

for  t  0,  e  e  E.  Note  that  the  conditional  distribution  (5.1)  is  inde¬ 
pendent  of  the  history  x  .  As  a  consequence,  it  is  clear  that  the  condi- 

n 

tional  distribution  G(dt;  x^,  e)  depends  on  x^  only  through  s^, 
l.e.,  G(dt;  x^,  e)  *  G(dt;  s^,  e).  (Recall  that  in  Section  4, 

G(dt;  x^,  e)  simplified  only  to  the  form  G(dt;  x^,  e).)  The  specific 
form  of  G(dt;  s,  e)  for  s  £  S,  e  e  E(8) ,  is  given  by  the  formula 

G(dt;  8,  e)  »  p(8,  e)  G(dt;  s) 


where 


G(dt;  s)  - 
q(a)  - 

p(8  ,e) 


qCs)  exp(-q(8)t)dt 


I  >^(e) 

eeECs) 


8,e 


X(e)  •  r  ^/qCs)  . 
s  ,e 


Hence,  it  follows  that 


(5.2) 


P{S_,  -  s'.  e  dt|X„  -  xj 


'irfl 


I  -  s’  e*  ,  -  e,  A__^,  e  dt,  X  -  x  } 

.  irt-1  '  n+1  ’  n+1  ’  n  n 

eeE(s  ) 
n 


'''S+l  ■  '•  Vl  '  ■  *n' 


I  p(8';  s  ,  e)  G(dt;  s  ,  e) 

eeE(8  )  " 

n 


p(8';  8^)  G(dt;  s^) 


where  p(8';  s)  ■  p(s' ;  s,  e)  pCs ,  e) .  It  easily  follows  from 

eeE(s) 

(5.2)  that  (S  ;  n  >  0)  is  a  Markov  chain  on  state  space  S  and  that  the 
n  — 

A  's  are  conditionally  independent  given  (S  :  n  >  0) ,  where 
n  n  — 


P{A  e  dtfS  :  m  >  0}  -  G(dt;  S  )  . 
irri  '  m  —  n 


A  well-known  characterization  of  continuous  time  Markov  chains  then  implies 
that  (S(t)  :  t  0)  is  the  (minimal)  Markov  process  corresponding  to  the 
generator  Q  ■  (Q(8p  s^)  :  82  e  S),  where 


25 


Q(8^,  s^)  ■  pCs^;  8^  q(8^)  ,  *1  ^  ®2 

Q(8j,  82)  “  “q(8j)  ,  «i  “  «2  * 

We  have  therefore  calculated  the  generator  of  the  contlnuous-tlae  Markov 
chain  associated  with  a  time- homogeneous  6SMP  in  which  all  clocks  are 
exponential. 

Perhaps  the  most  important  characteristic  of  a  continuous- time  Markov 
chain  Is  that  its  long-run  behavior  may  be  easily  calculated.  Specifi¬ 
cally,  if  the  internal  state  sequence  X  has  an  Invariant  probability 
distribution  n,  then  the  "Induced"  distribution  it(s)  ■  ii({s}  *  [0,®)  ^ 
C(a})  can  be  determined  as  the  probability  solution  of  the  system  of 
linear  equations  where  R(8,  s')  -  p(a’ ;  s) .  By  contrast, 

the  full  set  (2.12)  of  equations  for  n  typically  Involves  solving  an 
integral  equation. 

We  turn  now  to  semi-Markov  processes.  Consider  a  time-homogeneous 

GSMP  in  which  NCs';  s,  e)  -  ECs')  for  all  s',  s  e  S,  e  e  E(s).  In  this 

case,  for  n  >  1,  C  *0  for  e  e  E(S  ),  so  that 
—  n  ,e  n 

F<dt;  X  ,  e)  -  P  (dt)  a.s. 
n  e 

for  n  >  1.  As  a  result,  G(dt;  X  ,  e)  ■  G(dt;  S,,  e)  a.s.  for  n  >  1, 

—  n  n  — 

so  that  G(dt;  X  ,  e)  -  P{A  e  dt,  e*  ,  •  elx  }  depends  on  the 
o  xtTi  ttri  n 

history  of  X  only  through  S^.  Then,  for  n  1, 
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(5.3) 


«  I  G(dt;  S  ,  e)  p(8*;  s  ,  e) 

eeE(S  )  " 

n 

-  p(s’  ,  S  )  F(dt;  S  ,  s'  ) 
n  n 


where 


G(dt;  8,  e)  -  P{r“^^Y  e  dt}  n  p{r“^  ,Y  ,  >  t) 

e'eE(8)  ® 

e'l^e 


(Y  has  distribution  F  (•)) 
e  e 


pCs';  s)  -  J  G(®;  s,  e)  pis';  s,  e) 
ecEis) 


F(dt;  8,  s')  -  G(dt;  s,  e)  pCs';  s,  e)/p(8';8) 

-CE(8) 


It  follows  easily  from  (5.3)  that  (S  :  n  >  0)  is  a  Markov  chain  on  state 

n 

space  S  (having  transition  probabilities  p(s';  s))  and  that  the  A  's 

n 

are  conditionally  independent  r.v.'s  given  (S^  :  n  0) ,  where 

dt|S  :  m  >  0}  -  F(dt;  S  ,  S_^, )  . 
n+1  '  m  —  n  n+1 

By  definition,  (S(t)  ;  t  0)  is  therefore  a  semi-Markov  process.  As  in 
the  continuous-time  Markov  chain  context,  the  probability  it(8) 

“  7t({8}  [0,“)  ^  CCs))  can  be  calculated  as  the  solution  to  a  suitable 

system  of  linear  equations  (see  CINLAR  (1975)).  Thus,  the  analytical 
theory  for  the  steady-state  of  both  continuous- time  Markov  chains  and 
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serai-Markov  processes  is  considerably  simpler  than  that  encountered  for 
GSMP's. 

We  take  the  view  that  continuous- time  Markov  chains  and  semi-Markov 
processes  play  the  same  role  within  the  DEDS  area  that  linear  systems  play 
in  the  development  of  CVDS.  This  is  because  the  full  analytical  theory  of 
both  areas  can  only  be  brought  to  bear  in  the  specialized  settings 
mentioned  above. 

6.  REGENERATIVE  GSMP'S 

We  shall  now  show  that  an  important  class  of  time-homogeneous  GSMP's 
can  be  treated  as  regenerative  stochastic  processes.  Although  the  results 
to  be  described  here  are  far  from  the  most  general  possible,  they  are 
intended  to  give  a  flavor  of  what  can  be  expected  in  general. 

We  will  consider  GSMP's  in  which  the  distributions  F  (•)  have 

e 

certain  special  characteristics.  Suppose  that  a  distribution  F  satisfy¬ 
ing  F(0)  *  0  has  a  density  f  such  that  the  associated  hazard  rate 
h(t)  ■  f(t)/(l  -  F(t))  is  bounded  above  and  below  by  finite  positive 
constants  (i.e.,  inf{h(t)  :  t  0}  >  0,  8up{h(t)  ;  t  0}  <  ®).  The 
distribution  F  is  then  said  to  be  exponentially  bounded.  (We  use  this 
term  because  for  every  exponentially  bounded  distribution,  there  exist 
positive  constants  ^2  such  that  exp(-Xjt)  _<  1  ~  F(t)  ^  expC-X^t).) 

The  tall  distribution  function  of  a  positive  r.v.  can  be  represented 
in  terms  of  its  hazard  rate: 
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t 

1  -  F(t)  ■  exp(-  /  h(8)d8)  . 

0 

Hence,  if  F  is  an  exponentially  bounded  distribution  with  0  <  a  h(t) 

<  P  <  “, 

(6.1)  F(dt)  -  h(t)(l  -  F(t)dt 

a  exp(-Pt)dt 
6  •  p  exp(-Pt)dt 

where  6  -  a/p.  We  will  now  exploit  the  above  inequality  to  develop  a 
regenerative  structure  for  GSMP’s. 

The  idea  is  that  under  (6.1),  we  can  write 

(6.2)  F(dt)  =  6  •  P  exp(-pt)dt  +  (1~6)  0(dt) 

where  Q(dt)  is  a  probability  distribution  function.  Thus,  F(dt)  is,  with 
probability  6,  exponential  with  parameter  P.  Hence,  with  positive  prob¬ 
ability,  an  exponentially  bounded  clock  acts  like  a  memoryless  exponential 
clock.  This,  in  turn,  leads  to  regeneration. 

To  be  more  specific,  suppose  that  Fg(*)  is  exponentially  bounded 
for  all  e  e  E.  Let  a(e) ,  P(e)  be  the  lower  and  upper  bounds  on  the 
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associated  hazard  function  h^(*).  Then,  for  each  s  e  S,  we  have  the 
inequality 


(6.3) 


P{Sj  -  s' ,  e*  -  e,  c  dt (X^  -  (s,  t,  c) } 
>  £(8)p(8';  s,  e)  qCs)  exp(q(s)t)dt 


where 


£(8)  -  I  o(e)/p(e) 
e£E(8) 


pCs';  s,  e)  -  p(s';  s,  e)  P(e)r  /qCs) 

8 


q(s)  -  I  P(e)i 
eeE(s) 


s.e 


By  writing  the  inequality  (6.3)  as  an  equality  (in  the  same  fashion  as 
(6.2)  was  obtained  from  (6.1)),  we  see  that  if  Sq  ■  a,  then  with  probabil¬ 
ity  £(s),  (Sp  is  independent  of  Xq.  This  may  appear  to  suggest 

that  the  discrete-event  system  regenerates  with  probability  e(8)  every 
time  a  fixed  stated  s  £  S  is  hit.  Unfortunately,  this  reasoning  is 
invalid,  since  still  depend  on  C^.  Therefore,  we  need  to 

work  harder  to  obtain  regeneration. 

Assume  that  the  GSMP  satisfies: 

(6.4)  For  every  e  £  E,  there  exists  s  £  S  such  that  e  £  E(8)  and 

r  >0. 
s  ,e 

(6.5)  For  every  s,  s'  £  S,  there  exists  e,  s^  e^  ...,  9^,  e 
such  that 
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*  p(s';  8  ,  e  ) 
n  n 


>  0  . 


s  ,e 
n  n 


Condition  (6.5)  may  be  viewed  as  an  Irreduclblllty  hypothesis  on  the  GSMP. 
Under  assumptions  (6.4)  and  (6.5),  there  exists,  for  every  s,  s'  e  S, 


a  sequence  e ,  s  ,  e  ^  . . . ,  s^ ,  e^  such  that  the  GSMP  moves  from  s 
to  s'  with  positive  probability  through  the  Intermediate  states 

Sp  ...,  8^,  using  the  successive  trigger  events  . 

In  fact.  If  we  set  Sq,  ■  s,  8^^  -  s',  and  e^  ■  e,  we  have  the 
Inequality 


(6.6)  p{S^  *  ®i’  ®i  “  ®i-l  ’  ^1  ®  ‘**'l*  ^  ^  1^0  * 


>.  n  e(sj)  8^,  Cj^)  q(3^)  exp(-q(8j^)t^_j_^)dtj^^^ 

i«0 


In  fact,  conditions  (6.4)  and  (6.5)  allow  us  to  further  choose  the  path  so 
that  each  e  e  E(8q)  appears  In  the  set  {e^,  ...,  We  make  this 

choice  of  path  for  the  following  reason.  Note  that  given  Sq,  ..., 
e^,  ...,  the  clock  vector  Is  a  function  only  of  Cq, 

Ap  ...,  The  right-hand  side  of  (6.6)  shows  that  with  probability 

nj.g  e(8^)  pCs^j;  8^.  e^).  the  r.v.'s  Ap  ...,  may  be 

taken  as  Independent  exponential  r.v.'s.  Now,  with  our  choice  of  path,  we 
can  guarantee  that  for  each  e  e  ®  ^  ^(8^;  e^.j)  for 
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some  1  (1  <  1  <  nfl).  Then,  C.  *  0  and  it  follows  that  C  , ,  is  a 
function  purely  of  the  r.v.'s  •••»  observed  along  the  path. 

Since  the  are  independent  r.v.'s  with  positive  probability,  it 

follows  that  is  then  independent  of  (with  probability 

"i-o  “i-  'i''- 

For  each  s  e  S,  let  £<8)  -  c(s^)  e^l 

and  L(s)  be  the  length  of  the  path  constructed  above  (i.e.,  the  number  of 

states  visited).  We've  shown  that  if  S  “8,  then  with  probability  ^(s): 

m 

1)  (S  . , ,  A  . . ) .  1  <  1  <  L(8),  are  independent  r.v.'s  (S  is  the 

®t1  ittrl  —  —  ®T*1 

i'th  state  on  the  path  constructed  for  state  s,  and  A  ..  is  a 

fflTl 

corresponding  exponential  r.v.),  independent  of  X  , 

m 

2)  is  independent  of  X^. 

1)  and  2)  imply  that  ^^tttfL(s)+k  •  ^  2.  Independent  of  X^,  so  that 

the  sequence  ^^^nH-L(8)+k’  ^ittfL(8)+k^  :  k  0)  is  then  automatically 
independent  of  X^.  Combining  this  with  1),  we  conclude  that  with  proba¬ 
bility  e(S  ),  ( ( S ,  ,  A  )  :  k  >  1)  is  independent  of  X  .  Fix  a 
state  s  e  S.  We  have  just  shown  that  there  exists  a  random  subsequence 
T(l),  T(2),  ...  of  hitting  times  of  s  for  which  the  "cycles" 

{((Sj,  A^)  :  T(n)  <  j  _<  T(tH-l))  :  n  _>  1}  are  i.i.d.  The  random  subse¬ 
quence  of  regeneration  times  is  obtained  from  the  original  hitting  time 
sequence  by  flipping  a  coin  having  probability  of  success  eCs).  If  the 
coin  flip  is  successful,  then  the  next  LCs)  (S,  A)-tuples  are  generated 
using  the  algorithm  described  above.  This,  in  turn,  gives  rise  to  the 
desired  regenerative  structure.  We  have  thus  established  the  following 
result. 
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THBOtEM  4,  Consider  a  time- homogeneous  GSMP  satisfying  (6.4)  and  (6.5), 
for  which  Fg(*)  is  exponentially  bounded  for  all  e  e  E.  If  there 
exists  s  e  S  such  that  *  s  infinitely  often,  (S(t)  :  t  0)  is  a 
regenerative  process. 

An  interesting  feature  of  the  above  regenerative  construction  is  that 

while  the  r.v.'s  ((S  . .  ,  A  )  :  k  >  1)  are  Independent  of  X  ,  it  is 

nrrlc  nfrk  lo 

not  true  that  (X  :  k  >  1)  is  independent  of  X  .  Thus,  while  the 

ntTK.  —  111 

output  process  (S(t)  :  t  0)  is  regenerative,  the  internal  state 
sequence  may  not  be  regenerative.  A  similar  situation  arises  when  we 
consider  the  regenerative  structure  of  a  continuous-time  Markov  chain  from 
a  GSMP  viewpoint.  It  is  well  known  that  the  consecutive  times  at  which  the 
chain  hits  a  fixed  state  constitute  regeneration  times  for  the  associated 
(S,  A)  sequence.  On  the  other  hand,  the  full  vector  of  clock  read¬ 

ings  does  not  regenerate  at  such  hitting  times.  In  particular,  assuming 
all  speeds  are  unity,  the  differences  between  clock  readings  are  preserved 
from  one  transition  of  the  full  clock  sequence  to  the  next.  This  preserva¬ 
tion  of  memory  holds  even  at  transition  times  to  a  fixed  (physical)  state. 
Thus,  the  full  clock  vector  does  not  typically  regenerate,  even  when  the 
GSMP  is  a  continuous-time  Markov  chain. 

Suppose  |S(  <  ®  and  fix  s'  e  S.  From  (6.6),  it  follows  that 

(6.7)  P{T(s')  >  L  I  Xq  -  (s,  t,  c)}  >  e 

for  all  (s,  t,  c)  e  E,  where  e  -  min{e(s)  ;  s  £  S),  L  ■  fflax{L(s)  :  s  e  S), 
T(s')  ■  min{n  >  1  :  S  ■  s'}.  A  standard  "geometric  trials”  argument  then 
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proves  that  s'  is  visited  infinitely  often,  yielding  the  following 
corollary. 

C(M.OLLAKT  1.  Consider  a  time-homogeneous  GSMP  satisfying  (6. A)  and 
(6.5),  for  which  F^(*)  is  exponentially  bounded  for  all  e  e  E.  If 
(S(  <  ®,  then  (S(t)  :  t  0)  is  a  regenerative  process. 

A  regeneratalve  process  is,  in  some  sense,  a  stochastic  process  gener¬ 
alization  of  a  sequence  of  i.l.d.  r.v.'s.  As  a  result,  we  should  expect 
behavior  similar  to  that  typical  of  an  i.l.d.  sequence;  this  behavior 
Includes  strong  laws  and  central  limit  theorems. 

THEOREM  5.  Consider  a  time- homogeneous  GSMP  satisfying  (6.4)  and  (6.5), 
for  which  F^(*)  is  exponentially  bounded  for  all  e  e  E.  If  (s|  <  ®, 
there  exist  finite  (deterministic)  constants  r(f),  a(f)  such  that  for 
every  initial  distribution 

,  t 

—  /  f(S(u))du  r(f)  P  a.s. 

/  f(S(u))du  -  r(f))  ->  a(f)  N(0,1)  P  -  weakly 


as  t  ®. 


The  proof  of  Theorem  5  may  be  found  in  the  appendix.  An  important 
feature  of  Theorem  5  is  that  the  steady-state  limit  constants  r(f)  and 
(7(f)  are  Independent  of  p.  (Compare  with  Proposition  3.)  Note  also  that 
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if  we  view  f(S(t))  as  the  rate  at  which  cost  accrues  at  time  t,  then  the 
total  cost  C(t)  of  running  the  GSMP  up  to  time  t  has  a  distribution 
which  may  be  approximated  (in  distribution)  ’s 


C(t)  K  N(r(f)t,  cj^f)t)  . 

The  central  limit  theorem  of  Theorem  5  also  has  Important  implications  for 
steady-state  simulations  of  discrete-event  systems,  since  virtually  all 
steady-state  confidence  Interval  methodologies  (see  Chapter  3  of  BRATLEY, 
FOX,  and  SCHRAGE  (1987))  are  based  on  such  a  result. 

The  main  point  of  this  section  is  that  regenerative  ideas  can  be 
applied  to  discrete-event  systems.  The  construction  of  the  associated 
regeneration  times  is  typically  more  complicated  than  that  for  simpler 
processes,  such  as  continuous-time  Markov  chains.  Whenever  regenerative 
structure  is  present,  we  can  expect  results  similar  to  Theorem  5.  In 
addition  to  the  regenerative  structure  identified  here,  HAAS  and  SHEDLER 
(1987a) ,( 1987b)  have  identified  regeneration  (of  a  different  character)  in 
a  number  of  other  GSMP  contexts.  Thus,  we  view  the  laws  of  large  numbers 
and  central  limit  theorems  described  here  as  being  typical  of  a  large  class 
of  discrete-event  systems. 

7.  LIKELIHOOD  KATIOS  FOE  GSMP'S 

Let  A  S  and  suppose  that  we  wish  to  calculate  p{S(A)  <  t},  where 

S(A)  ■  lnf{t  0  ;  S(t)  c  A).  Typically,  this  probability  needs  to  be 

numerically  calculated;  simulation  Is  generally  the  most  popular  numerical 
approach. 


35 


In  many  situations,  we  expect  that  P{S(A)  ^  t}  Is  small.  For 
example,  If  A  Is  the  set  of  "failed  states"  of  a  discrete-event  reliabil¬ 
ity  system,  then  p{s(A)  ^  t}  will  be  small  if  the  system  is  reliable. 
Unfortunately,  naive  simulation  is  highly  inefficient  for  such  problems; 
many  replications  will  be  necessary  in  order  for  the  system  to  experience  a 
reasonable  number  of  failures. 

A  powerful  technique  that  can  be  used  in  such  situations  is  importance 
sampling.  The  idea  is  to  simulate  the  system  so  as  to  bias  it  toward 
failure;  the  estimator  must  then  be  altered  so  as  to  compensate  for  the 
"biased  dynamics."  The  adjustment  factor  needed  is  called  a  likelihood 
ratio. 

Consider  a  GSMP  of  the  type  described  in  Section  3.  The  probability 

distributions  that  govern  the  dynamics  of  the  system  are  the  conditional 

distributions  p( • ;  x  ,  e)  and  P(*;  x  ,  e).  Let  P(*)  denote  the 

n  n 

probability  distribution  of  the  internal  state  sequence  X  under  these 
conditional  distributions,  and  let  E(»)  be  the  corresponding  expectation 
operator. 

To  perform  importance  sampling,  we  need  to  specify  the  alternative 
conditional  distributions  that  will  appropriately  bias  the  dynamics  of  the 
system.  For  x^,  e,  let  JgC*;  e)  and  F(»;  x^,  e)  be  condi¬ 
tional  distributions  having  the  property  that  there  exist  functions 
q(  • ;  X  ,  e)  and  f  ( • ;  x  ,  e)  such  that 
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(7.1) 


p(*;  e)  »  q(*;  x^,  e)  x^,  e) 

(7.2)  F(dt;  X  ,  e)  ■  f(»;  x  ,  e)  F(dt;  x  ,  e)  . 

n  n  ~  n 

denote  the  probability  distribution  and  expectation 
operator  corresponding  to  the  conditional  distributions  p(*;  x  ,  e)  and 
F( • ;  x^,  e) .  The  following  result  is  a  straightforward  generalization 
of  the  likelihood  ratio  ideas  in  GLYNN  and  IGLEHART  (1988b). 

THEOREM  6.  Consider  a  GSMP  with  conditional  distributions  p( * ;  x^,  e) , 
£(•;  e)  ,  F(dt;  x^,  e) ,  F(dt;  x^,  e)  satisfying  (7.1)-(7.2). 

Let  T  be  a  stopping  time  relative  to  the  internal  state  sequence  X 
(i.e.,  1(T  -  n)  is  a  function  of  X^) ,  and  let  Y  *  ^(^0*  **’ 
real-valued.  Then, 


E  YI(T  <  <»)  •  E  YI(T  <  »)L_ 
p  T 


(the  equality  should  be  interpreted  as:  if  one  expectation  exists,  then 
both  do  and  they  are  equal) ,  where 


L^  -  g(Aj;  X^_j,  e*)  q(S^;  X^_j ,  e*) 

and 


g(t;  X  e)  -  f(r  t;  x  ,  e)  n  F(r  ,t;  x  ,  e')/F^(r  ,t;  x  ,e'), 

n  8,e  n  e'eE(8)  " 

e'  j*e 
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Theorem  6  is  the  key  to  importance  sampling  for  GSMP’s.  Rather  than 

replicate  copies  of  the  r.v.  YI(T  <  ®)  under  to  estimate 

a  »  E  Y1(T  <  “),  we  can  replicate  copies  of  Y1(T  <  “)L_  under  P  to 
p  T  ~|i 

estimate  a.  By  choosing  P^  appropriately,  significant  Improvements  in 
computational  efficiency  over  conventional  simulation  can  be  achieved. 

Likelihood  ratio  ideas  can  also  be  applied  to  parameter  optimization 
of  discrete-event  systems.  Specifically,  likelihood  ratio  methods  can  be 
used  to  obtain  an  efficient  means  of  estimating  the  gradient  of  the  objec¬ 
tive  function  via  simulation  (see  GLYNN  (1987)).  This,  in  turn,  can  be 
used  to  develop  a  simulation-oriented  gradient-based  algorithm  for  optimiz¬ 
ing  discrete-event  processes. 

The  likelihood  ratio  methods  described  here  are  but  two  examples  of 
how  the  GSMP  structure  of  a  discrete-event  system  can  be  used  to  obtain 
computational  enhancements  to  numerical  algorithms  for  discrete-event 
systems. 

8.  APPENDIX 

Proof  of  Theorem  2.  We  first  show  that  X  is  weakly  continuous  on  the 

state  space  I,  i.e.,  if  x  ,  x  e  Z  and  x  x  as  n  ®,  then  P(x  ,  •) 

■>  P(x,  •)  as  n  ®,  where  ■>  denotes  weak  convergence.  (Recall  that 

P(x,  •)  •  PtX  ,  e  •  (  X  *  x}  is  the  transition  function  of  X.) 
n+l  n 

Let  (V(e,  Cg)  :  e  e  E)  be  a  collection  of  independent  r.v.'s  having 
marginal  distributions  specified  by 
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P{V(e,  c  )  >  t}  -  F  (t  +  c  )/F  (c  )  . 
e  e  e  e  e 

Then,  under  (4.1)  and  (4.2),  the  conditional  probability  distribution 
G(*;  x^,  e)  =*  G(*;  x^,  e) ,  where 

G(u;  (s,t,c),  e)  «  Ef(V(e* ,  c^,)  :  e'  e  E(8)) 

and  f(v(e')  :  e'  e  E(s))  *  1(  min  v(e')/r  ,  <  u, 

e'eE(s)  s.®  - 

argmln  v(e')/r  ,  *  e). 

e'eE(s).  8.® 

By  conditions  11)  and  ill).  It  is  evident  that  F  (t+c)/F  (c)  is 

e  e 

continuous  In  c  at  every  t.  Hence,  (V(e,  c*)  :  e  e  E) 

e 

*>  (V(e,  c  )  :  e  e  E)  whenever  (c*  :  e  e  E)  (c  :  e  e  E).  Since  f 

e  e  e 

is  continuous  at  (V(e,  c^)  :  e  e  E)  (we  use  il)-iv)  here).  It  follows 

that  G(*;  x' ,  e)  *>  G(*;  x,  e)  whenever  x'  x.  Thus,  the  distributions 

of  the  trigger  event  e*  and  A  are  (weakly)  continuous  In  x  ■  (s,  t,  c). 

Consequently,  the  distribution  of  Xj^  »  (Sj^,  A,  Cj^)  is  (weakly)  continuous 

in  Xq  =  (s,  t,  c) ,  thereby  proving  the  required  continuity  of  X. 

The  second  step  involves  showing  that  for  every  e  >  0,  there  exists  a 

compact  subset  F  of  2  such  that  P(x,  F)  1  -  e  uniformly  in  x  e  Z. 

Let  r  »  min{r  :  s  e  S,  e  e  E(s)}.  Select  K,  using  condition  v) ,  so 
s  ,e 

that  Fg(i'K  +  c)/Fg(c)  <  e/|Ej  uniformly  in  e  and  c.  Put 
F  ■  {x  -  (s,  tc)  e  Z  :  t,  £  K,  e  e  e}  and  observe  that 
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P(x,  <  I  p{v(e,  c  )/r  ^  >  K} 

eeE(s)  ®  ®’® 

<  I  F  (rK  +  c)/F  (c)  <  e  . 
eeE(s)  ®  ® 

Let  p”(x,  •)  be  the  n-step  transition  function  of  X  (l.e.  , 
P^(x,  •)  “  ^  *  f  ^0  *  ®o“®  X  e  2,  set 


n-1 


I  P^(x,  •)  . 


j“0 


Observe  that  ^  ^  ^  tight. 

Prohorov's  theorem  asserts  the  existence  of  a  subsequence  n*  and  a 
probability  it  on  Z  such  that  ■>  x,  A  standard  argument  (see,  for 

example,  KARR  (1975))  then  uses  the  weak  continuity  of  X  to  prove  that  x 
Is,  In  fact,  Invariant  for  X. 


Proof  of  Theorem  3.  We  first  note  that  since  |s|  <  °°,  f  is  bounded  and 
it  suffices  to  prove  only  that  ^  *•  ®y  the  statlonarity  of  X  under 

X ,  this  is  equivalent  to  showing  that  ^-^^2  ^  **  will  follow  if  we 

can  prove  that  E{A2(Xj  ■  x)  is  uniformly  bounded  In  x. 

For  s,  s' ,  e  e  E(8),  let 


M(8' ;  s,  e) 
m(8';  8,  e) 


min  V(e',  0)/r  ,  ,  , 

e'eN(8*;8,e)  ’ 

EM(8';  8,  e)  . 
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G.early,  mCs';  s,  e)  <  J  t  F  (dt)/r,  where  r  Is  the  minimum  over 

~  eeE  [0,®)  ® 

the  speeds  r  ,  ,  of  condition  ill).  To  finish  the  proof,  note  that 

s  ,e 

Cj  g  *■  0  for  e  e  N(S^;  Sq,  e^)  and  hence  E{A2|Xj  ■  x)  ^  Em(S^;  Sq,  e^) 

<  I  /  t  F  (dt)/r. 

eeE  [0,“) 

Proof  of  Theorem  5.  We  need  to  verify  the  hypotheses  of  the  regenerative 
strong  law  of  large  numbers  and  central  limit  theorem  (see  SMITH  (1955)). 
Fix  s'  e  S,  and  let  t  be  the  first  m>^L(s')  such  that  the  next  L(s') 
states  to  be  visited  (after  m)  form  the  specified  path  for  s' .  (We  will 
base  our  regenerations  on  visits  to  s'.) 

Since  |s|  <  “,  f  Is  bounded  so  that  It  suffices  to  verify  that  the 
moment 


e{(  I  -  x} 

1*1 

is  bounded  In  x.  We  first  note  that  there  exists  L  such  that 

sup  P{t  >  L|Xq  ■  x}  <1 

X 

(use  (6.7)).  Since  the  tail  of  x  is  then  geometrically  dominated,  this 
Implies  that  E{x^|Xq  ■  x]  is  bounded  In  x,  for  all  P  >  0.  Now 
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e{(  I  A)2(Xq-x} 

1,1  ^ 

£  e{x^  max  jXj.  ■  x} 

1  <  1  <  T  ^  ^ 

<  -  x}  E^^^{  max  aJ(X«  -  x} 

^  I  <  i  < 

<  E^^^{x^|x  -  x)  E^^^{  I  aJ(Xq  -  x}  , 

"  1-1  ^  ^ 

so  we  need  only  show  the  second  factor  above  Is  bounded  In  x.  But 

e{  I  aJIXq  -  x}  -  i  E{A^  I(T  >  1) |Xq  -  x} 

1-1  ^  1-1 

CO 

<  I  E^^^{A®|X  -  x}  •  >  i|x  -  x)  . 

1-1  ^  " 

Again,  since  the  tall  of  t  Is  geometrically  dominated,  It  suffices  to 

O 

prove  that  E{A^(X^  -  x}  Is  bounded  in  x.  But  for  x  -  (s,  t,  c) ,  we 
can  select  e  e  E(8)  so  that  r  >0,  Then, 
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00 


Since 

which 


EUfiXj, 


“  x}  <  J 


0  ‘  -  W 


4-1/® 

c  +r  t 

®  e  8,e 

-  /  expl-  /  h^(u)du)dt 

0  c. 


1  /ft 

<  /  exp(-a(e)  r  t  )dt  . 
-  0 


(s|  <  *,  it  is  evident  that  •  x)  is  bounded  in  x,  from 

it  follows  that  E{A®(Xq  -  x)  is  bounded  in  x. 
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